/**
 * @file   effective_ensemble.hpp
 * @author Kimmo Luoma <kimmo@PHOTON3>
 * @date   Fri Sep  7 14:16:36 2012
 * 
 * @brief  Effective ensemble of NMQJ method for a single time instant. 
 *Contains the relevant states
 *for the system dynamics and their probabilities and defines some associated
 *methods. 
 * 
 * 
 */
#ifndef EFFECTIVE_ENSEMBLE_HPP
#define EFFECTIVE_ENSEMBLE_HPP
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <complex>
#include <map>
#include <utility>
using namespace boost::numeric::ublas;
class effective_ensemble{
private:
  std::map<int, vector<std::complex<double> > > states;
  int ntot;
protected:
  std::map<int,int> occ;
public:
  effective_ensemble(){};
/** 
 * copy constructor 
 * 
 * @param eff reference to object type effective_ensebmle
 */
  effective_ensemble(const effective_ensemble& eff);
  
  /** 
   * constructor 
   * 
   * @param mat d x n matrix<std::complex<double> >  of the states.
   * d is the dimenstion of Hilbert space and n is the number of the 
   * states. 
   * @param vec vector<double> of probabilities of length n.
   * @param o 
   *
   **/
  effective_ensemble(matrix<std::complex<double> > mat,vector<int> o);
  /** 
   * Copy constructor
   * 
   * @param eff object of type effective ensemble.
   * 
   * @return reference to object type effective ensemble.
   */
  effective_ensemble& operator=(const effective_ensemble eff);
  /** 
   * Obtain a state vector in the effective ensemble.
   * 
   * @param alpha index of the state vector
   * 
   * @return vector<std::complex<double> a state vector 
   */
  vector<std::complex<double> > get_state(int alpha);
  /** 
   * Set a particular state vector.
   * 
   * @param alpha index of state vector
   * @param vec state vector <std::complex<double> >
   */
  void set_state(int alpha,vector<std::complex<double> > vec);
  /** 
   * Return a probability of a state in the ensemble.
   * 
   * @param alpha index of the particular state.
   * 
   * @return probability of \f$\vert \psi_\alpha \rangle\f$ in the ensemble.
   */
  double P(int alpha);
  /** 
   * Return the number of states in the ensemble.
   * 
   * 
   * @return number of states in the ensemble type int.
   */
  int get_dim();
  /** 
   * Return a probability vector.
   * 
   * 
   * @return vector<double> containing probabilities of the different states
   * in the ensemble.
   */
  vector<double> P();
  /** 
   * Return the density matrix of the system 
   * \f$ \rho=\sum_\alpha p_\alpha\vert\psi_\alpha\rangle\langle\psi_
   * \alpha\vert\f$.
   * @return matrix<std::complex<double> > rho
   */
  matrix<std::complex<double> > rho();
  /** 
   * Return the number of states in the ensemble.
   * 
   * 
   * @return number of states in the ensemble type int.
   */
  int get_number_of_states();
  /** 
   * Return density matrix as a vector \f$\vec{\rho}\f$.
   * mapping is \f$\vec{\rho}_{i*d+j}=\rho_{i,j}\f$, where
   * \f$d\f$ is the dimension of the Hilbert space.
   * 
   * @return vector<std::complex<double> >
   */
  vector<std::complex<double> > rho_vec();
  /** 
   * Get occupation of state alpha
   * 
   * @param alpha index of the state
   * 
   * @return <int> number of ensemble members in state \f$\vert \psi_\alpha \rangle\f$.
   */
  int occupancy(int alpha);
  
  /** 
   * Calculate the variance 
   * 
   * 
   * @return matrix of variances
   */
  vector<double> variance();

  /** 
   * Set occupancy of state alpha
   * 
   * @param alpha 
   */
  void set_occupancy(int alpha,int n);
  
  /** 
   * Ge all occupation numbers
   * 
   * 
   * @return vector<int> 
   */
  vector<int> occupancy();
  
};

#endif
